use ../processed/NLSY79.dta, replace

gen AFQT3_a = max(AFQT3-1/3,0)
gen AFQT3_b = max(AFQT3-2/3,0)
gen female = (sex==2)
gen age_a = max(age-35,0)
gen age_b = max(age-45,0)
gen finc_a = max(finc_pct-1/3,0)
gen finc_b = max(finc_pct-2/3,0)

local ctrl locE locU urban age age_a age_b nsibs

// parental education
forvalues j=1/5{
	gen hgc_m`j' = (hgc_m==`j')
	gen hgc_f`j' = (hgc_f==`j')
	local ctrl `ctrl' hgc_m`j' hgc_f`j'
}
// cohort
local yob_all
forvalues jj=58/64{
	gen yob`jj' = (yob==19`jj')
	local yob_all `yob_all' yob`jj'
}
// Census division
forvalues j=2/9{
	gen divisiona`j' = (divisiona==`j')
	local ctrl `ctrl' divisiona`j'
}

local xlist
forvalues j=9/18{
	gen x`j' = (hgc_r>=`j')
	local xlist `xlist' x`j'
}

gen xc = max(hgc_r-12,0)
local ctrl0 AFQT3 AFQT3_a AFQT3_b `ctrl'

local ctrl1 `ctrl0' `yob_all' female blk hsp // cross-sectional sample
local ctrl2 `ctrl0' yob62 yob63 yob64 female blk hsp finc_pct finc_a finc_b
local ctrl3 `ctrl0' `yob_all' blk hsp // feale only
local ctrl4 `ctrl0' `yob_all' blk hsp // male only
local ctrl5 `ctrl0' `yob_all' female hsp // black or hispanic
local ctrl6 `ctrl0' `yob_all' female // non-black, non-hispanic
local ctrl7 `ctrl0' `yob_all' female blk hsp // direct cost only
local ctrl8 `ctrl0' `yob_all' female blk hsp // indirect cost only

local iv pub4 tui4 locE17 locU17 // direct cost only
forvalues j=1/6{
	local iv`j' `iv'
}
local iv7 pub4 tui4 
local iv8 locE17 locU17

local cond1 `"if cs_smpl==1"'
local cond2 `"if !missing(finc_pct)"'
local cond3 `"if female==1"'
local cond4 `"if female==0"'
local cond5 `"if hsp==1 | blk==1"'
local cond6 `"if hsp==0 & blk==0"'
local cond7
local cond8

matrix A = J(8,21,.)
mat colnames A = ols ols_se iv iv_se gap gap_se cw cw_se tw tw_se me me_se fs_pub4 fs_pub4_se fs_tui4 fs_tui4_se fs_locE fs_locE_se fs_locU fs_locU_se N
forvalues j=1/8{
	qui reg hgc_r `iv`j'' `ctrl`j'' [pw=wgt] `cond`j''
	predict z, xb
	if ( `j'<=6 | `j'==7 ){
		matrix A[`j',13] =  _b[pub4]
		matrix A[`j',14] = _se[pub4]	
		matrix A[`j',15] =  _b[tui4]
		matrix A[`j',16] = _se[tui4]
	}
	if ( `j'<=6 | `j'==8 ){
		matrix A[`j',17] =  _b[locE17]
		matrix A[`j',18] = _se[locE17]	
		matrix A[`j',19] =  _b[locU17]
		matrix A[`j',20] = _se[locU17]	
	}
	matrix A[`j',21] = e(N)
	
	qui reg hgc_r `ctrl`j'' [pw=wgt] `cond`j''
	predict x_res, resid
	qui foreach v in z xc `xlist'{
		regress `v' `ctrl`j'' [pw=wgt] `cond`j''
		predict `v'_res, resid
	}
	local ctrlx
	local ctrlxc
	foreach v in `ctrl`j''{
		gen x_`v' = hgc_r*`v'
		gen xc_`v' = xc*`v'
		local ctrlx `ctrlx' x_`v'
		local ctrlxc `ctrlxc' xc_`v'
	}
	local xlist_r
	foreach v in `xlist'{
		local xlist_r `xlist_r' `v'_res
	}
	
	// predict E[Y|X,W]
	regress logwage hgc_r `ctrlx' `ctrl`j'' [pw=wgt] `cond`j''
	predict vhat1, resid
	gen mhat1 = _b[hgc_r]*x_res
	foreach v in `ctrl`j''{
		replace mhat1 = mhat1 + `v'*_b[x_`v']*x_res
	}
	regress z_res hgc_r `ctrlx' `ctrl`j'' [pw=wgt] `cond`j''
	predict zhat1, xb
	regress logwage `xlist'  `ctrlxc' `ctrlx' `ctrl`j'' [pw=wgt] `cond`j''
	predict vhat2, resid
	gen mhat2 = 0
	foreach v in `ctrl`j''{
		replace mhat2 = mhat2 + `v'*( _b[x_`v']*x_res+_b[xc_`v']*xc_res ) 
	}
	foreach v in `xlist'{
		replace mhat2 = mhat2 + _b[`v']*`v'_res
	}
	regress z_res `xlist'  `ctrlxc' `ctrlx' `ctrl`j'' [pw=wgt] `cond`j''
	predict zhat2, xb	
	
	corr x_res z_res [aw=wgt] `cond`j'', covar
	scalar m_xx = r(Var_1)
	scalar m_xz = r(cov_12)
	
	// OLS
	reg logwage hgc_r `ctrl`j'' [pw=wgt] `cond`j'', vce(cluster cty_grp)
	scalar b_ols = _b[hgc_r]
	matrix A[`j',1] = b_ols
	matrix A[`j',2] = _se[hgc_r]
	predict e_ols, resid
	
	// IV
	ivregress 2sls logwage (hgc_r=z) `ctrl`j'' [pw=wgt] `cond`j'', vce(cluster cty_grp)
	scalar b_iv = _b[hgc_r]
	matrix A[`j',3] = b_iv
	matrix A[`j',4] = _se[hgc_r]
	predict e_iv, resid
	
	// IV-OLS difference 
	gen temp = e_iv*z_res/m_xz - e_ols*x_res/m_xx
	reg temp [pw=wgt] `cond`j'', vce(cluster cty_grp)
	drop temp
	matrix A[`j',5] = b_iv - b_ols
	matrix A[`j',6] = _se[_cons] 
	
	// covariate weight difference
	ivregress 2sls mhat1 (hgc_r=z) `ctrl`j'' [pw=wgt] `cond`j'', vce(cluster cty_grp)
	scalar b_c = _b[hgc_r]
	predict e_c, resid
	gen temp = ( e_c*z_res + vhat1*zhat1 )/m_xz - e_ols*x_res/m_xx
	reg temp [pw=wgt] `cond`j'', vce(cluster cty_grp)
	drop temp
	matrix A[`j',7] = b_c - b_ols
	matrix A[`j',8] = _se[_cons]  
	
	// treatment-level weight difference
	ivregress 2sls mhat2 (hgc_r=z) `ctrl`j'' [pw=wgt] `cond`j'', vce(cluster cty_grp)
	scalar b_ct = _b[hgc_r]
	predict e_ct, resid
	gen temp = ( (e_ct*z_res + vhat2*zhat2) - (e_c*z_res + vhat1*zhat1) )/m_xz
	reg temp [pw=wgt] `cond`j'', vce(cluster cty_grp)
	drop temp
	matrix A[`j',9] = b_ct - b_c
	matrix A[`j',10] = _se[_cons]  
	
	// marginal effect difference
	gen temp = ( e_iv*z_res - (e_ct*z_res + vhat2*zhat2) )/m_xz
	reg temp [pw=wgt] `cond`j'', vce(cluster cty_grp)
	drop temp
	matrix A[`j',11] = b_iv - b_ct
	matrix A[`j',12] = _se[_cons]  	
	
	drop z x_res z_res xc_res `xlist_r' `ctrlx' `ctrlxc' mhat1 zhat1 vhat1 mhat2 zhat2 mhat2 vhat2 e_ols e_iv e_c e_ct
}

clear
svmat A, names(col)
gen model = ""
replace model = "Cross-sectional sample only" in 1
replace model = "Family income sample" in 2
replace model = "Females" in 3
replace model = "Males" in 4
replace model = "Black or Hispanic" in 5
replace model = "Non-black, non-Hispanic" in 6
replace model = "Direct cost only" in 7
replace model = "Indirect cost only" in 8
order model, first
export delimited using ../result/decom_alt.csv, replace

